Pregunta 1:

Explique cual es la diferencia fundamental entre la estadística bayesiana y frecuentista.

La estadistica frecentista tiene una interpretación de la probabilidad como la frecuencia con que un evento ocurre en el tiempo, es decir, se requiere de repetir un experimento hasta converger a una probabilidad por grandes números. La estadistica bayesiana tiene una interpretación sobre la probabilidad que tiene más que ver con el conocimiento que se tiene sobre un determinado evento, es decir, se combina conocimiento posterior y datos obtenidos.

Ambas estadsticas tienen como base su interpretación de la probabilidad, lo cual hace que esta sea la diferencia fundamental entre estas.

Pregunta 2:

Discuta la siguiente afirmación La inferencia bayesiana permite fácilmente utilizar distintos tipos de información.

La probabilidad bayesiana, toma en cuenta distintos tipos de información a la hora de calcular la distribución de probabilidad de un parámetro. Principalmente se toma en cuenta la probabilidad previa (puede ser conocimiento experto o una hipotesis) y la probabilidad que se obtiene de los datos. En conjunto, ambos tipos de información sirven para obtenet una probabilidad de salida, la cual puede ser vista como una atualización de la probabilidad previa.

La probabilidad previa (prior probability) puede tener una distribución distinta a la de los datos, esto facilita el manejo de hipotesis.

Pregunta 3:

Explique la diferencia entre prior probability y posterior probability

prior probability: Corresponde a la probabilidad que se cree o se tiene antes de evaluar los datos. Esta puede ser una suposición o tambien conocimiento experto. Esta probabilidad forma parte de la inferencia para obtener la posterior probability.

posterior probability: Corresponde a la probabilidad que se obtiene luego de evaluar los datos. Esta es la probabilidad que se busca obtener.

En resumen, su diferencia radica en como se obtienen y la función que cumplen en la inferencia bayesiana.

Pregunta 4:

El estadista bayesiano “Bruno Finetti” menciona la siguiente frase en su libro de probabilidad: La probabilidad no existe. Lo que en verdad quizo decir es que la probabilidad es un método para describir incertidumbre en un observador con conocimiento limitado. Discuta esta información utilizando el ejemplo del lanzamiento del globo terraqueo visto en clases. ¿Que significa decir “la probabilidad de que sea agua es 0.7”?

La frase “la probabilidad de que sea agua es 0.7” no significa que este sea un valor constante que se encontró y siempre será así. La probabilidad tiene una distribución que depende del conocimiento previo y de los datos obtenidos.

En consecuencia, esta frase significa que hasta ahora, en base al conocimiento previo y a los datos, se puede estimar este valor de probabilidad, el cual viene dado por el conocimiento actual y está asociado a certeza de ocurrencia más que un valor determinado.

Pregunta 5:

¿ Que ventaja entrega que la distribución de la posterior este en la misma familia de distribución de probabilidad que la del prior?. De un ejemplo de alguna distribución que posee este comportamiento.

Si ambos, prior y posterior, tienen igual distribución, eso quiere decir que en cada iteración se irá ajustando la curva. Gracias a esto se verá una convergencia de los parametros de la distribución. En un caso contrario, donde prior y posterior tienen distribuciones distintas, la convergencia se complica y el resultado no necesariamente será de una distribución conocida.

Por ejemplo, si se tiene una distribución gaussiana para prior probability y posterior probability, lo que ocurrirá es que los parámetros \(\sigma\) y \(\mu\) serán los que se ajustarán.

Pregunta 6:

Señale y explique los pasos de la grid approximation para obtener la posterior y responda las siguientes preguntas:

1. ¿Cual de los pasos señalados nos permite obtener una distribución de la posterior mas precisa?.
2. Cuales son las limitaciones de la grid approximation.

A continuación se explican los pasos de la grid approximation

• 1. Definir grilla: Definir la discretización a realizar
• 2. Evaluar grilla: Calcular el prior para cada valor de la grilla
• 3. Calcular likelihood en para cada grilla
• 4. Multiplicar cada prior con su likelihood (para cada valor de grilla)
• 5. Dividir el arreglo resultante en la suma total, es decir, normalizar

1. En el primer paso se decide la discretización, es decir, la precisión de la grilla con respecto a la curva real.

2. Cuando se tiene un gran número de variables este método se vuelve inviable, ya que se debe hacer grilla por cada parámetro del sistema. Esto tiene un limite computacional por debajo de otros métodos.

Pregunta 7:

¿ Por qué es necesario aprender a trabajar con muestras de la posterior?.

Debido a que el método más utilizado y con mejor desempeño para calcular aproximaciones de varios parámetros, multiniveles y más complejas de la Posterior es el método de Monte Carlo usando Cadenas de Markov. Y esto es importante ya que éste método no devuelve la distribución de posterior directamente, si no que precisamente devuelve una muestra.

Pregunta 8:

Señale si las siguientes afirmaciones son verdaderas o falsas, en el caso que sean falsas justifique su respuesta:

• Los point estimates de la posterior no entregan información relevante en un estudio.

• Un intervalo de confianza es un intervalo dentro del cual un valor de parámetro no observado cae con una probabilidad particular, mientras que un intervalo de credibilidad es un rango de valores en el que se estima que estará cierto valor desconocido.

• La principal ventaja de HPDI frente a un intervalo de credibilidad es que si la posterior no distribuye de forma normal, el HPDI será capaz de detectar los puntos de interés, mientras que un intervalo de credibilidad lo ignoraría al asumir simetría.

• [F] Los point estimates de la posterior no entregan información relevante en un estudio.

Falso, si bien es cierto que los point estimates reducen un montón la información relevante de un estudio, dado que proporcionan información bastante comprimida del fenómeno que se está estudiando, los point estimates nos pueden ofrecer información importante de manera resumida antes de inspeccionar con mayor ahínco los datos. El máximo estimado de la posteriores, por ejemplo, nos daría un primer acercamiento a valores importantes del modelo, como la simetría de la distribución. Pese a lo anterior, siempre es recomendable dar la mayor cantidad de información posible de la distribución, y no quedarse sólo con los point estimates, si no también con otros estimadores, como funciones de densidad, HPDI, etc..

• [F] Un intervalo de confianza es un intervalo dentro del cual un valor de parámetro no observado cae con una probabilidad particular, mientras que un intervalo de credibilidad es un rango de valores en el que se estima que estará cierto valor desconocido.

Con respecto al intervalo de confianza falta agregar que este intervalo aparece luego de repetir infinitamente la sample distribution del expertimento, y este cae con cierta frecuencia dentro del rango. Para el intervalo de credibilidad, no es sobre cierto valor desconocido, sino más bien una cantidad específica de la probabilidad de Posterior.

• [V] La principal ventaja de HPDI frente a un intervalo de credibilidad es que si la posterior no distribuye de forma normal, el HPDI será capaz de detectar los puntos de interés, mientras que un intervalo de credibilidad lo ignoraría al asumir simetría.

Pregunta 9

Suponga que tiene dos especies de pandas. Cada una de las especies es igual de común y es imposible distinguirlas físicamente. Una de las diferencias entre las especies es el tamaño de sus familias. Si denotamos por \(\theta\) a la especie del panda se tiene que, cuando la especie es \(\theta = 1\) tiene dos bebes un \(10\%\) de las veces mientras que la especie \(\theta = 2\) tiene dos bebes un \(20\%\) de las veces, mientras que el resto de veces ambas especies tienen un solo bebe.

Suponga que usted esta intentando determinar la especie de un panda que que tiene como registro de nacimientos al conjunto \(D\), considere quela especie de un panda que acaba de dar a luz a dos bebes, es decir \(D = \{\text{dos bebes}\}\)

• ¿Cual es la probabilidad de que pertenezca a la especie \(1\)?

• Suponga ahora que el mismo panda acaba de dar a luz y esta vez es solo un bebe. Calcule la probabilidad posterior de que el panda sea de especie \(1\). ¿Que cambia con el calculo anterior?

• Suponga que le ofrecen hacer un test genético a su panda, como suele ser común con los test no es perfecto y le mencionan las siguientes características:

  • La probabilidad de que correctamente idenfitique a la especie \(1\) es de \(0.8\)
  • La probabilidad de que correctamente identifique a la especie \(2\) es de \(0.65\)

Se administra el test y se obtiene un resultado positivo a la especie \(1\). Sin utilizar la información en \(D\) calcule la probabilidad posterior de que su panda sea de la especie \(1\). Repita sus cálculos utilizando la información recopilada en \(D\). ¿En que varían sus resultados?

Usando el teorema de Bayes, nos enfocamos en obtener los parámetros del teorema:

“Cada una de las especies es igual de común”:

\[\mathbb{P}(\theta_1) = \mathbb{P}(\theta_2) = 0.5\] Se calcula la probabilidad de cada tipo de panda dado los datos. \[\mathbb{P}(D|\theta_1) = 0.1,\, \mathbb{P}(D|\theta_2) = 0.2\]

Se calcula la probabilidad de los datos, sin información adicional. \[ \mathbb{P}(D) = \mathbb{P}(D|\theta_1) \cdot \mathbb{P}(\theta_1) + \mathbb{P}(D|\theta_2) \cdot \mathbb{P}(\theta_2) = 0.5 \cdot 0.1 + 0.5\cdot0.2 = 0.15\] Por ello, usando el teoréma Bayesiano:

\[\mathbb{P}(\theta_1|D) = \frac{\mathbb{P}(D|\theta_1) \cdot \mathbb{P}(\theta_1)}{\mathbb{P}(D) } = \frac{0.5 \cdot 0.1}{0.15} = 0.\overline{3}\]

Usando la misma metodología anterior:

\[ \mathbb{P}(D|\theta_1) = 0.1 \cdot 0.9 = 0.09, \, \mathbb{P}(D|\theta_2) = 0.2 \cdot 0.8 = 0.16\] \[\mathbb{P}(D) = \mathbb{P}(D|\theta_1) \cdot \mathbb{P}(\theta_1) + \mathbb{P}(D|\theta_2) \cdot \mathbb{P}(\theta_2) = 0.5 \cdot 0.09 + 0.5\cdot0.16 = 0.125\] \[\mathbb{P}(\theta_1|D) = \frac{\mathbb{P}(D|\theta_1) \cdot \mathbb{P}(\theta_1)}{\mathbb{P}(D) } = \frac{0.5\cdot 0.09}{0.125} = 0.36\]

Los resultados muestran que aumenta la probabilidad de que sea un panda tipo \(\theta_1\), dado que este nuevo trozo de información lo favorece.

Es posible que el test haya obtenido un resultado positivo de dos maneras: - Correctamente eligió a \(\theta_1.\, \mathbb{P} \to 0.8\). - Incorrectamente no eligió a \(\theta_2.\, \mathbb{P} \to (1 - 0.65) = 0.35\). Y viceversa para el \(\theta_2\).

Siguiendo el mismo procedimiento bayesiando de las preguntas anteriores:

\[\mathbb{P}(\theta_1) = 0.5, \,\mathbb{P}(\theta_2) = 0.5\]

\[\mathbb{P}(T|\theta_1) = 0.8, \,\mathbb{P}(T|\theta_2) = 0.65\] \[\mathbb{P}(T) = \mathbb{P}(T|\theta_1) \cdot \mathbb{P}(\theta_1) + \mathbb{P}(T|\theta_2) \cdot \mathbb{P}(\theta_2) = 0.5 \cdot 0.8 + 0.5\cdot 0.65 = 0.725\]

Entonces la probabilidad de que sea \(\theta_1\) es:

\[\mathbb{P}(\theta_1|T) = \frac{\mathbb{P}(T|\theta_1) \cdot \mathbb{P}(\theta_1)}{\mathbb{P}(T) } = \frac{0.8 \cdot 0.5}{0.725} = 0.551..\] Ahora se repetirán los cálculos considerando la información anterior, esto significa utilizar el Posterior como prior en nuestros cálculos:

\[\mathbb{P}(\theta_1) = 0.\overline{3}, \,\mathbb{P}(\theta_2) = 0.\overline{6}\]

\[\mathbb{P}(D|\theta_1) = 0.1, \,\mathbb{P}(D|\theta_2) = 0.2\]

\[\mathbb{P}(D) = \mathbb{P}(D|\theta_1) \cdot \mathbb{P}(\theta_1) + \mathbb{P}(D|\theta_2) \cdot \mathbb{P}(\theta_2) = 0.8 \cdot 0.\overline{3} + 0.65 \cdot 0.\overline{6} = 0.7\]

\[\mathbb{P}(\theta_1|D) = \frac{\mathbb{P}(D|\theta_1) \cdot \mathbb{P}(\theta_1)}{\mathbb{P}(D) } = \frac{0.8 \cdot 0.\overline{3}}{0.7} = 0.38\] Calculando el Posterior de \(\theta_1\) usando \(D\), sabiendo previamente los resultados del test, que favorecen a \(\theta_1\), el resultado muestra que \(\theta_1\) es ahora más favorable que antes.

Preguntas: Introducción a Grid Approximation

Las primeras dos preguntas de esta tarea tienen como objetivo introducirlos en la inferencia Bayesiana utilizando la técnica Grid Approximation para obtener una aproximación de la posterior. Al finalizar los problemas ustedes deberán ser capaces de visualizar los efectos que tiene el prior en la posterior, saber cómo realizar una Grid Approximation y comprender como utilizar Percentile Interval (PI) en una posterior.

Pregunta 1:

Considere el dataset “moneda.csv” donde se encuentran los resultados de un experimento lanzando una moneda, el objetivo de esta pregunta es estudiar mediante técnicas de inferencia Bayesiana el valor de la probabilidad de que salga cara, representado por el valor \(1\). Puede usar la librerira rethinking durante toda esta pregunta (si lo desea).

dataMoneda <- read.csv("moneda.csv", header = TRUE)

• Programe el metodo grid approximation para distintos tamaños de experimento. ¿Como van cambiando las curvas posterior?

• Repita el mismo análisis anterior pero utilizando el método de Laplace (no necesita programar el método, puede utilizar la libreria “rethinking”). ¿Como se comparan con los resultados anteriores?.

• Grafique la densidad de la posterior y encuentre la proporción de los siguientes defined boundaries:

  • \((0, 0.4)\)
  • \((0.4,0.7)\)
  • \((0.7,1)\)

¿Como puede interpretar los resultados?

• Calcule un intervalo de credibilidad al \(50\%\), \(75\%\) y \(95\%\). ¿Como se puede interpretar los resultados? ¿Cual podría ser un problema al usar intervalos de credibilidad?.
• Genere los intervalos HDPI para \(50\%\), \(75\%\) y \(95\%\), compárelos con los intervalos de credibilidad de la parte anterior. ¿En que se diferencian los intervalos de credibilidad con los HDPI?.

grid_posterior <- function(wins, total, grid_size){
  
  grid <- seq(0, 1, length.out = grid_size)
  prior <- rep(1, grid_size)
  likehood <- dbinom(wins, total, prob=grid)
  posterior <- likehood * prior
  posterior <- posterior / sum(posterior)
  #print(posterior)
  list(posterior = posterior, grilla = grid)
}


total_list <- c(5, 10, 20, 50, 100, 200, 500, 1000)
win_list <- c()

for (total in total_list){
  wins <- length(dataMoneda[1:total,][dataMoneda[1:total,] == 1])
  win_list <- append(win_list, wins)
}

grid_detail = 900
Posterior <- c()
index_list <- c(1:length(total_list))
for (i in index_list){
  res <- grid_posterior(win_list[i], total_list[i],  grid_detail)
  print(plot(res$grilla, res$posterior, type="b"))
  
}

## NULL

## NULL

## NULL

## NULL

## NULL

## NULL

## NULL

## NULL

A medida que más datos se van agregando a los cálculos, se observa que la distribución del Posterior se va consolidando en una opción más plausible y por lo tanto la curva está más concentrada en un valor más denso.

library(rethinking)


total_list <- c(5, 10, 20, 50, 100, 200, 500, 1000)
win_list <- c()

for (total in total_list){
  wins <- length(dataMoneda[1:total,][dataMoneda[1:total,] == 1])
  win_list <- append(win_list, wins)
}

grid_detail = 900
Posterior <- c()
index_list <- c(1:length(total_list))
for (i in index_list){
  globe.qa <- quap(
  alist(
    W ~ dbinom(W+L, p), 
    p ~ dunif(0, 1)
    ), 
  data=list(W=win_list[i], L=(total_list[i] - win_list[i])))
  precis(globe.qa)
  sample.quap <- extract.samples(globe.qa)
  #hist(sample.quap)
  dens(sample.quap)
  #res <- grid_posterior(win_list[i], total_list[i],  grid_detail)
  #print(plot(res$grilla, res$posterior, type="b"))
}

Con el método de Laplace, observamos el mismo comportamiento, si bien la librería rethinking muestra los resultados de una manera que más difícil de visualizar a lo largo de las iteraciones (la densidad se empieza a agrupar en intervalos cada vez más pequeños, de la misma manera que con el método Grid approximation)

globe.qa <- quap(
  alist(
    W ~ dbinom(W+L, p), 
    p ~ dunif(0, 1)
    ), 
  data=list(W=win_list[8], L=(total_list[8] - win_list[8])))
  precis(globe.qa)

  sample.quap <- extract.samples(globe.qa)
  #hist(sample.quap)
  dens(sample.quap)

  prop_00_to_04 <- sum(sample.quap < 0.4)/length(sample.quap[, 1])
  prop_04_to_07 <- sum(sample.quap > 0.4 & sample.quap < 0.7)/length(sample.quap[, 1])
  prop_07_to_10 <- sum(sample.quap > 0.7 & sample.quap < 1.0)/length(sample.quap[, 1])
  
  #prop_07_to_1 <- c(quantile(sample.quap[, 1], 0.7), quantile(sample.quap[, 1], 1.0))

prop_00_to_04

## [1] 0

prop_04_to_07

## [1] 1

prop_07_to_10

## [1] 0

Estos resultados expresan directamente que el 100% de la probabilidad del posterior (de que la moneda salga cara) se encuentra entre 0.4 y 0.7

print("Intervalo de credibilidad al 50%")

## [1] "Intervalo de credibilidad al 50%"

PI(sample.quap[, 1], prob=0.5)

##       25%       75% 
## 0.5411294 0.5624577

print("Intervalo de credibilidad al 75%")

## [1] "Intervalo de credibilidad al 75%"

PI(sample.quap[, 1], prob=0.75)

##       12%       88% 
## 0.5334151 0.5699446

print("Intervalo de credibilidad al 95%")

## [1] "Intervalo de credibilidad al 95%"

PI(sample.quap[, 1], prob=0.95)

##        3%       98% 
## 0.5208478 0.5830639

La información que nos devuelve en este caso el intervalo de credibilidad, es que la mayoría de la densidad se encuentra concentrada entre el 0.52 y 0.58, lo cual es un rango de probabilidad bastante pequeño. Esto igual tiene sentido dado que se ha tirado una moneda (que lo más probable es que esté ligeramente cargada) tantas veces que existe bastante certidumbre sobre el Posterior más probable.

El problema de los intervalos de credibilidad es que no se comportan bien cuando la probabiliad del posterior es asimétrico, esto porque puede que no contengan el valor más plausible dentro de sus rangos. Afortunadamente en este dataset, el comportamiento de la moneda está modelado de manera bastante simetrica respecto a las dos salidas, por lo que no vemos este problema.

print("HPDI al 50%")

## [1] "HPDI al 50%"

HPDI(sample.quap[, 1], prob=0.5)

##      |0.5      0.5| 
## 0.5403592 0.5616133

print("HPDI al 75%")

## [1] "HPDI al 75%"

HPDI(sample.quap[, 1], prob=0.75)

##     |0.75     0.75| 
## 0.5333602 0.5697956

print("HPDI al 95%")

## [1] "HPDI al 95%"

HPDI(sample.quap[, 1], prob=0.95)

##     |0.95     0.95| 
## 0.5198888 0.5817645

Los valores obtenidos con este método más estable frente a densidades no simétricas son bastante parecidos a los anteriores, lo que nos demuestra de manera más analítica que la densidad de probabilidad del Posterior es simétrica o tiene un nivel muy bajo de \(skewness\).

df = read.csv("no+mordidas.csv")
head(df)

df_2 <- c(df[, 1], df[, 2])
mordidos <- df_2[df_2 == 1]
cant_mordidos <- length(mordidos)
cant_no_mordidos <- length(df_2) - cant_mordidos
cant_muestra <- length(df_2)

p2_res <- grid_posterior(cant_mordidos, cant_muestra, 10000)
grilla <- p2_res$grilla
post   <- p2_res$posterior
index <- which.max(post)
max_posterior <- grilla[index]
#print(plot(p2_res$grilla, p2_res$posterior, type="b"))

print(paste("Maximo valor de la posterior:" , as.character(max_posterior)))

## [1] "Maximo valor de la posterior: 0.496049604960496"

pred_list = c(1:10000)
for (i in c(1:10000)){
  
  pred_list[i] <-  rbinom(1, size=500, prob=max_posterior)
}
hist(pred_list)

# Calculamos cuantos fueron mordidos en el primer mes dado que fueron mordidos en el segundo, y su inverso.
primer_mes_dado_mordido_segundo <- df[, 1][df[, 2] == 1]
mordidos_mes1_dado_mordido_mes2 <- length(primer_mes_dado_mordido_segundo[primer_mes_dado_mordido_segundo == 1])
no_mordidos_mes1_dado_mordido_mes2 <- length(df[, 1]) - mordidos_mes1_dado_mordido_mes2

# Calculamos el posterior usan grid approximation
p3_res <- grid_posterior(mordidos_mes1_dado_mordido_mes2, length(df[, 1]), 10000)
p3_grilla <- p3_res$grilla
p3_post   <- p3_res$posterior
p3_index <- which.max(p3_post)
p3_max_posterior <- p3_grilla[p3_index]

# Generamos la simulacion en base a ese posterior 10000 veces
p3_pred_list = c(1:10000)
for (i in c(1:10000)){
  
  p3_pred_list[i] <-  rbinom(1, size=length(df[, 1]), prob=p3_max_posterior)
}
hist(p3_pred_list)

# Leer información
data_notas <- read.csv("notas.csv")

# Función para crear likelihood dado mu y sigma
grid_function <- function(mu,sigma){
   return(sum(dnorm(data_notas$Notas, mean=mu, sd=sigma)))
}

# Valores de la grilla
grid_mu <- seq( from=5 , to=7, length.out=100 )
grid_sigma <- seq( from=0.5 , to=1.5 , length.out=100 )

# Se crea la grilla 2d
data_grid <- expand_grid(grid_mu,grid_sigma)

# Se guarda la likelihod
data_grid$likelihood <- map2(data_grid$grid_mu,data_grid$grid_sigma, grid_function)

# Se transforma el forma de map2 a una columna
data_grid <- unnest(data_grid,cols = c("likelihood"))

# Valores de los priors
prior_mu <- rep(1/length(grid_mu), times = 100)
prior_sigma <- rep(c(2/(3*length(grid_sigma)/2), 1/(3*length(grid_sigma)/2)), each = length(grid_sigma)/2)

# Se crea la grilla 2d de priors
prior <- expand_grid(prior_mu,prior_sigma)

# Se calculan los valores del prior
data_grid$prior <-  map2(prior$prior_mu,prior$prior_sigma, prod)
data_grid <- unnest(data_grid,cols = c("prior"))

# Se calcula el posterior
data_grid$unstd_posterior <-  data_grid$likelihood * data_grid$prior

# Se estandariza el posterior
data_grid$posterior <- data_grid$unstd_posterior/sum(data_grid$unstd_posterior)

# Se ajustan los valores de la posterior para que no sean valores tan pequñeos
data_grid$posterior <- (data_grid$posterior - min(data_grid$posterior))/(max(data_grid$posterior)-min(data_grid$posterior))

# Punto de referencia
# Se recomienda cambiar estos valores por unos adecuados que le permitan estudiar
# Los valores de la distribución de mejor manera
valor_x <- 6.1
valor_y <- 1

# Grafico

punto_comparacion <- tibble(x = valor_x, y = valor_y)

plt <- data_grid %>%
  ggplot(aes(x = grid_mu, y = grid_sigma)) +
  geom_raster(aes(fill = posterior),
    interpolate = T
  )+ 
  geom_point(x = valor_x, y = valor_y, size = 1.3,color="white")+
  geom_label(
    data = punto_comparacion, aes(x, y),
    label = "Punto Comparación",
    fill = "green",
    color = "black",
    nudge_y = 0.1, # Este parametro desplaza la caja por el eje y
    nudge_x = -0.2 # Este parametro desplaza la caja por el eje x
  )+
  scale_fill_viridis_c() +
  labs(
    title = "Posterior para Mean y Standard Deviation",
    x = expression(mu ["Mean"]),
    y = expression(sigma ["Standar Deviation"])
  ) +
  theme(panel.grid = element_blank())

plt

# Codificamos los datos
x <- 1:length(data_grid$posterior)

# Sampleamos los indices
posterior_samples_aux <- sample(x,size = 1e4, replace = T, prob = data_grid$posterior)

# Obtenemos los verdaderos valores de la sampling distribution
posterior_samples <- data_grid[posterior_samples_aux,]

# Obtenemos solos los valores relevantes para la densidad
df <- data.frame(posterior_samples$grid_mu,posterior_samples$grid_sigma)

# Realizamos las densidades
dens(df)